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ABSTRACT 



Accordling to the theory of Kozai resonance, the initial mutual inclination 
between a small body and a massive planet in an outer circular orbit is as high 
as ~ 39.2° for pumping the eccentricity of the inner small body. Here we show 
that, with the presence of a residual gas disk outside two planetary orbits, the 
inclination can be reduced as low as a few degrees. The presence of disk changes 
the nodal precession rates and directions of the planet orbits. At the place where 
the two planets achieve the same nodal processing rate, vertical secular resonance 
would occur so that mutual inclination of the two planets will be excited, which 
might trigger the Kozai resonance between the two planets further. However, in 
order to pump an inner Jupiter-like planet, the conditions required for the disk 
and the outer planet are relatively strict. We develop a set of evolution equations, 
which can fit the N-body simulation quite well but be integrated within a much 
shorter time. By scanning the parameter spaces using the evolution equations, 
we find that, a massive planet (10M/) at 30AU with 6° inclined to a massive 
disk (50Mj) can finally enter the Kozai resonance with an inner Jupiter around 
the snowline. And a 20° inclination of the outer planet is required for flipping 
the inner one to a retrograde orbit. In multiple planet systems, the mechanism 
can happen between two nonadjacent planets, or inspire a chain reaction among 
more than two planets. This mechanism could be the source of the observed 
giant planets in moderate eccentric and inclined orbits, or hot-Jupiters in close- 
in, retrograde orbits after tidal damping. 

Subject headings: Celestial mechanics planetary systems: protoplanetary disks 
planets and satellites: dynamical evolution and stability-planets and satellites: 
formation 
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Introduction 



Kozai mechanism is a kind o: 



systems (ILidov 



1962 



Kozai 



1962 



secular effect tha t occurred in hierarchical three-body 



Naozetal 



2011bJ). In the limit of circular restricted 



three-body model, the test particle in an inner orbit can be pumped to a highly eccentric 



or inclinec 
is > 39.2° 



orbit as long as its initial inclina tions relative to the outer massiy e 



Kozai 



1962 



Innanen et al. 



pertu rber 



1997T ) . Further more, iLithwick fc Naozl (120111 ) show 
that, when the massive perturber is in an eccentric orbit, the effect of octupole terms in 
the perturbing function will be effective so that the inner test particle may be repelled to 
retrograde orbits relative to the massive planet orbit. 

One of the most prominent applications for Kozai mechanism is on the formation of 



orbital confi gurations 



(RM) effect (IRossiter 



or h o t Jupiters (HJ). R ecent observations of Rossiter-McLaughlin 



1924 



McLaughli 



19241 ) show that most of the HJs might be in 



orbits misaligned with stellar spins. Actually, fo r the 53 HJs with RM effect measurem ents 



at least 8 H 


s mi 


Brown et al. 


2012 



s might be in retrograde mot ions (IWinn et al 



Albrecht et al. 



2010 



Triaud et al. 



2010 



20121 ). As the classical core accretion scenario says that 
planets were formed in a protoplanetary disk surrounding the protostar, the existence 
of HJs in highly inclined orbits infers that some dynamical mechanisms to pump their 
i nclinations must be exis t after their for mation. As the so-called disk migration scenario 



( )Lin fc Papaloizou 



1986 



Lin et al. 



1996|) failed to explain the existenc e of HJs in retrograd e 



orbits, Kozai mechnism was invoked to excite the orbital inclination s (jWu fc Murray 



Fabrycky fc Tremaine 



2007 



Naoz et al 



2011a 



Nagasawa et al 



2008|). 



2003 : 



Wu fc Murravi fl2003h : 



Fabrycky fc Tremaind (120071 ) proposed that a third massive body 



(either a binary or a brown dwarf companion) with a high orbital inclination(> 39.2°) can 
trigger the Kozai resonance so that the orbital eccentricity of inner planets can be pumped 
up to near 1, which can be damped at periastron of the orbit, with its excited inclination 



-4- 



being preserved. However, the population studies establis h that on 



explained by Kozai migration due to binary companions (IWu et al. 



10% of HJs can be 



20071 ). while most of 



the HJ systems did not find any stellar or substellar companions. An alternative choice is 
whether the outer perturber can be replaced by a massive planet. Altho ugh this is possible 



a very high mutual inclination between the two planets is required. E.g. 



Naoz et al. 



(2011a) 



presents a flipping example with a 3Mj planet as t he outer perturber , while the initial 



Lithwick fc Naozl (j201lf ) shows that 



mutual inclination of the two planets is up to 71.5°. 
if the outer perturbers are in more eccentric orbits, the relative inclination can be reduced, 
but still as high as ~ 60° for retrograde motion to be occurred. Thus, the origin of such 
high mutual inclination itself merits explanations. 

In this paper, we propose a mechanism to efficiently excite planetary eccentricities and 
inclinations with an outer residual gas disk. After gas giants have formed and swept away 
the inner part of gas disk, the residual gas disk outside will perturb the architecture of inner 
planet systems. Due to the gravity of the residual disk, vertical secular resonances would 
occur between the very massive outer planet and the inner ones at some certain locations. 
Then the mutual inclination between the planetary orbits would be pumped. At this time, 
if the outer planet has a non-zero inclination relative to disk midplane, which might result 
from the previous planetary scattering, the mutual inclination is possible to raise up to the 
Kozai critical value, then the Kozai effect between planets would be induced. As a result, 
the eccentricities and inclinations of the inner planets would be excited to very high values. 



The effect o 



Nagasawa et al. 



gas disk in exciting 



(120031 ) 



Terquem et al. 



planet a ry eccentricities was als o studied by 



(120 id ) 



Teyssandier et al. 



( 120121 ) . etc. Our present 



work focus on how, with the aid of a disk, a two planet system will execute secular 
resonances between them in order to trigger the subsequent Kozai effect. We will also 
present the parameter studies with a set of evolution equations. The paper is organized as 
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follows: in section 2, we introduce the model and two examples, and the pumping region 
is displayed by scanning the a^o — ^2,0 plane. Then we point out the pumping mechanism 
is the secular resonances coupled with the following Kozai resonance, and calculate the 
location of secular resonances in the situation of small eccentricities and inclinations by 
timescale comparisons in section 3. In section 4, we deduce the changing rates of some 
crucial parameters in pumping process relative to any inertial plane, and compare them with 
N-body simulation results. In section 5, influence of planetary parameters are investigated. 
According to that, we give the critical pumping conditions for a fixed gas disk. Section 
6 displays situations in systems with more than two planets. Finally, discussions and 
conclusions are presented in section 7. 



2. Model and Examples 



We consider a planet system with two giant planets (denote as m,\ and m 2 for inner 
and outer planet, respectively) orbiting around a centr al star, with a protoplanetary disk 



whose inner part had been swept out by giant planets (IZhang et al. 



2008 



2010a 



Zhang fc Zhou 



bf). The gas disk is assumed to be a two-dimensional circular annulus for simplicity, 
with its mass distributed on the midplane. As both the mass and the angular momentum 
of the disk are much larger than those of the planets, we further suppose the gravity of 
the planets has no influence on the disk, i.e., the disk is invariable. As the disk exerts the 
gravity onto the planets, the equation of planet motion can be written as follows 

N 

f 1 ll) . 



d 2 n 



dt 1 



G(m + m i )( r 
r 



) + E Gm 



r, - n 



•..13 



11 



where is the position vector of the planets relative to the star, and 

d(j) 



$ = -G / £(r)rdr / 

JR in JO 



r 2 + r 2 — 2rr p cos sin Op) 1 / 7 



(2) 
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displays the gravitational potential from the disk ( iTerquem et al.ll2010l ). Ri n ,R ont is the 
inner and outer border of the disk. (a p , <p p , 9 P ) is the spherical coordinates of a planet in the 
coordinate system settled by the star and the disk midplane. S(r) is the mass density of 
the disk, and we use the most commonly exponential density distribution of the disk radius 
r, E(r) = £o(r/-R ou t)~ a . Total mass of the disk is settled by Mdisk and the expression for 
S is shown in Appendix B. 

We apply Runge-Kutta-Fehlberg 7(8) integrator to integrate Equations ([T]). Figure 
[T] gives a typical example, whose initial conditions are listed in table 1. We set the star 
mass m = 1M Q . The inner and outer boundary of the out gas disk are taken arbitrarily 
within the scope of disk observations. In order to satisfy the assumption that the angular 
momentum of the disk is overwhelming, we set the mass of the g as disk as 0.05M ffi . Though 



it is much larger than the average mass (O.O1M ) estimated by 



Williams fc Ciezal (1201 lh . 



it is still within th e reasonable range acc ording to the recent transitional disk observation 



(such as LkCa 15 (IKraus fc Ireland 



20121 )). The mass of outer planet is moderately bigger 



than the inner one for facilitating the excitation procedure, and the particular influence 
will be discussed in section [51 We take the initial eccentricity and inclination of the inner 
planet very small just to show the pumping mechanism. Eccentricity of the outer planet is 
set very small in order to conveniently compare with the results of the evolution equations 
(section HJ), and the non-zero eccentricity situation will be discussed in Section [5j 

We can see from Figure [1] that the inclination of the inner planet relative to disk 
midplane (Ii) goes up to near 50° within 0.3Myr. After around 0.4Myr, the eccentricity of 
the inner planet e\ begins to rise, accompanied with the mutual inclination between planets 
(/tot) declining. Ascending nodes of the two planets precess with same rates for most of 
the first 0.3Myr, which implies that it is the secular resonance that raises l\ and hence 
J tot . This triggers the whole excitation procedure. Argument of pericenter of the inner 
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planet u\ keeps librating during the cause. We also notice that for a while after 0.7Myr, I\ 
becomes larger than 90° and meanwhile e\ is close to 1. It provides a good opportunity for 
the planet to turn into a retrograde hot- Jupiter after considering tidal damping due to the 
central star. As comparison, the case with the same initial conditions except for free of gas 
is shown in the right. Eccentricities from planetary secular perturbations merely are much 
smaller, and mutual inclination keeps around 30° all the time. 

Figure [2] gives another example with smaller initial inclination of m 2 {J%$ = 10°) (The 
subscript means the initial value, and hereafter). And the mutual inclination of two 
planets could also be stirred up to 40° companied by the approaching nodes precession rates 
of the two planets. Then e± is pumped by Kozai effect with the sign of u\ librating. 

The initial inclination J 2 ,o an d semi-major axes are critical parameters for pumping 
occurring, so we scanned the phase space of a 10 — i2,o> an d for every case, extracted the 
maximum of 7 to t, h, ei and e 2 (short by J t ot,max, A,max, ei imax and e 2 , m ax hereafter) during 
the evolutions within lMyr. Figure [3] in filled color shows those with parameters the same 
as Figured] except the variable a 10 and 7 2 ,o- Both Ji imax and I to t,max have obvious minimums 
at around a 10 = 3.5au when / 2 = 0°. And the area above the contour line of Itot,max = 40° 
coincides with that above the line of ei jmax = 0.1 (the discrepancy in their upper left corner 
is due to a longer Kozai timescale than lMyr, so there is no enough time for e\ to rise), 
which implies that eccentricity pumping is attributed into Kozai effect after inclinations 
have been excited. In Figure [3jd, the eccentricity of m 2 is also excited in the region that 
the planetary secular resonances occur (when J 2 ,o < 35°). Although e 2jmax becomes large in 



some regions either due t o the secular resonan c e or combined wit 



i the Kozai oscillations 



(-^2,0 > 35°) from gas disk tjTerquem et al. 



2010 



Teyssandier et ajj 



20121 ) . it still maintains 



less than 0.1 in most cases, which is the basis of simplifications in the derivation of the 
evolution equations in section HI 
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These pumping cases represent a possible scenario to excite efficiently the eccentricities 
and inclinations of planets when planets are far away from each other. And the pumping 
critical angle is much lower than the Kozai critical angle because of the initial inclination 
excitation. From the nearly equal rates of change of nodes of two planets we have deduced 
it is secular resonance that excites the inclinations. And we will further verify that in the 
next section by frequency and timescale comparisons. 



3. Conditions for secular resonances (ESR and VSR) 



Secular evolution dominates dynamics of a planetary system when planets are far 
away from the star and they are not close to any low-order mean-motion resonances. In 
this context, once the precessio n frequencies of plane t s are integer multiples of each other, 



secular resonance would occur (ILithwick &: Wu 



2011 



Nagasawa fc Idal |2000| ) .At the place 



where the timescales of perihelion (nodal) processing rate of two planets are equal due to 
the disk and mutual planetary perturbations, secular resonance would occur, which are 
called as eccentric (vertical, respectively) secular resonance, and shorted as ESR (VSR, 
receptively). 

In order to obtain the timescales more explicitly, we first assume the initial eccentricities 
and inclinations of both planets are small before they are effectively excited. We further 
assume that m 2 > m 1 , and a 2 ^> a 1; so the evolution of m 2 is dominated by perturbations 
from the disk, and that of m x is mainly affected by perturbations from m 2 (also see Figure 

ED. 



Under these assumptions, we use Lagrange equations (Murray fc Dermott 



19991) to 



derive the apsidal and nodal precession rate exerted by the disk gravity (see Appendix B 
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for details) 

3 

O i)disk = — Kcos/j, (3) 
2 

Wj.disk = K, (4) 

rii 

where Jj and Qi are the inclinations and ascending nodes of the two planets (i = 1,2) with 
respect to the disk midplane, Ui is the argument of perihelion, rij is the angular velocity of 
planetary mean motion, and 

-a + 2 -l + ^GMd^ 
1-7T-+ 2 -1-a 2it* ut W 

is merely related to the disk parameters, 77 = Ri Q /R out , a is the exponential index of disk 
profile. 

In deriving Equations (J3j)-(SJ) , terms with e 2 and sin 2 / have been eliminated to simplify 
the expressions, which is suitable before the exciting of e and /. Meanwhile, under these 
assumptions, the semi-major axis a, eccentricity e and inclination I of each planet have no 
secular trend from disk gravity (see Appendix B). So the timescales of planet apsidal and 
nodal precession from disk gravity can be estimated by 2ir/u and 2n/Cl separately. Then 
the timescales of the outer planet are 

2vr 2tt 
T n 2 = 7- , r W2 = . (6) 

"2,disk W 2 ,disk 



19991 ) to obtain 



Moreover, we apply the secular perturbation theory (Murray fc Dermottl 
the precession timescale of the inner planet due to planetary interactions. There are two 
eigenfrequencies gi, #2 (where g\ > g 2 ) for e — u solution and one eigenfrequency / for I — Q 
solution in two-planet systems. So 

27T 27T 

= ~r, V = — (<) 

/ 9i 

can be used to display the precession of Qi and u\ respectively. 
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Figure H] shows these timescales versus the inner planet's semi-major axis, with initial 
condition the same as Figure [3] except 7 2 ,o = 0. Tn 1 and tq 2 , t Ui and r U2 respectively have 
one cross point in Figure SJ The x-coordinations of the points display the value of ai o when 
VSR and ESR occur, which roughly match the location of pumping at ^2,0 = in Figure [3j 
And the y-coordinations estimate the timescales of sec ular resonances, which are much less 



than the average ages of gas disk (IHaishch et al. 



2001b- 



When J 2 ,o > 0, tq 2 becomes larger(Equation (jSJ), then the cross point of m 1 and r^ 2 
would move inward along the tq 1 line. So it only provides the estimation of the inner 
border of the excitation region in Figure [21 In order to estimate the excitation region more 
precisely, we will give the evolution equations of the elements in the next section. 



4. Evolution equations at arbitrary inclinations 



To obtain the quantitative description of planetary orbits when secular resonance 
happens, we will develop a set of simplified equations to describe the evolutions of 
ex,Ux,Ii, Hi, 12,^2, which are suitable for arbitrary inclinations (but still require for small 
e 2 ). We set the disk midplane as the reference plane, which is assumed to coincide with the 
equatorial plane of the center star. So our derivations are different with iNaoz et al.l (l2011bl ) 
in the context of three-body systems, as their reference plane is the invariable plane of the 
system. 



At first, according to 



Mardling fc Lin 



(120021 ). the secular evolution of the elements of 



rax effected by m 2 is expressed into the angular momentum vector h= rxr, the Runge-Lenz 
vector e and q = hxe (see Equation (1A3[) - (1A6h ) (The hat indicates the unit vector). And 
for mi-, the correspond vectors are H, E and Q. Then time-averaging is executed, first over 
the inner orbit for eliminating eccentric anomaly Ex then over the outer orbit for removing 



- 11 - 



E 2 (the results see Equation flA9j) - flA14l) ). Then after, we expand the two groups of unit 
vectors (e, q, h) and (E, Q, H) into terms with J l5 ui, Qi and J 2 , u 2 , fi 2 separately (see 
Equation ( 1A15P ). This is the key step to make the final formulas relative to an arbitrary 
plane rather than the invariable plane of two orbits. Finally, we obtain the evolutions of the 
elements due to planetary perturbation up to the quadrupole terms, without any reductions 
on the eccentricities and inclinations (see Equation (IA16j) - (lA21|) ) . It is worth mentioning 
that, the evolutions of e\ and u\ has no assum ption of AQ = 7r, so has more terms than the 



quadrupole parts in formula (C9) and (C5) of 



Naozetal 



f|2011bh . 



The disturbing from gas disk is been considered independently. Details are in Appendix 
B. The final evolutions can be acquired by adding the two parts together, 



doc , doc >, / doc 



x represents the six elements Ii,i2,fii,fi 2 ,ei an d We set e 2 = as e 2 keeps small in most 
cases (Fig. EJi), then the six equations presented by the above one become closed (we call 
them "the evolution equations" hereafter). 

We made comparisons for the two parts of the evolution equations by drawing 
\og[(dx/dt)p/(dx/dt)<iisk\ from true N-body simulation in Figure [5] As was expected, for the 
elements of mi, (dx/dt) p ^> (dx/dt)disk in most time, and for f2 2 , (<if2 2 /<it)p (rf^/c^disk 
all the time. As for J 2 , the influence from disk is much smaller because of the small e 2 (see 
Equation ( 1B8I) ). These can be utilized in the further deductions and simplifications. 

Via the evolution equations, we can calculate the evolution of inclination and 
eccentricity of the inner planet more quickly. The dashed line in Figure [1] displays the 
integration results from the evolution equations. Except for some delay, both the trend and 
the amplitude are fitted pretty well. And further, we use the evolution equations to make 
scanning over a 10 — / 2 ,o (the black contour lines in Figure |3k-c, which is made up of the 
maximums of the evolutions of lMyr or before e\ > 0.99 for every case) to compare with 



the full N-body results. The simplified results agree well qualitatively with the full N-body 
ones, except some malposition, which mostly results from the quadrupole approximation 
for disk gravity. 



5. Parameter analysis 

According to Figure [HJ Kozai resonance between planets occurred above the contour line 
of -^tot,max = 40°, and the retrograde motion of mi happened above the line of Ii imax = 90°. 
As planets were thought to be coplanar at their earliest stage, lower values of extremum of 
these two contour lines would make the generation of retrograde motion easier. For this 
purpose, we investigate the dependence of the minimums of h tm3X = 90° and itot,max = 40° 
on G^o and m 2 with the evolution equation (jSJ) (Fig. [6]), with the parameters the same as 
Figure H] except for the variables a 2) o, rn 2 and the scanned 1^,0, ai,o- Filled color contour is 
composed of the values of y-coordinate of the extremum (/2,o,min) ; which means the smallest 
inclination of rri2 for the onset of Kozai effect (itot.max = 40°) or for m\ retrograding 
(A,max = 90°). The solid line contour is built up by x-coordinates of the extremum, which 
signify the locations of m\ when VSR between planets will occur. 



Considering a . 
for a IMq star, see 



upiter-mass pla net most probably formed outside the snowline (2.7au 



Ida k, Lin 



20041 ). we constrain the interesting scope beyond 2.7au for 
solid contour in Figure El We can see that with a O.O5M gas disk ranging from 50au to 
lOOOau, a Jupiter-mass planet at ~ 2.7au will be pumped by Kozai effect with a 5mj planet 
at 25au and inclined > 10° relative to the disk midplane. Further more, it can be flipped 
into a retrograde orbit by a giant planet at ~ 25 au with mass of 5m j and inclination 
> 30°, or mass of 10m j with inclination > 20°. As every a^o — I2.0 scanning involved in 
Figure [6] is made up of the < 1 Myr integrations of the evolution equations, a general gas 
disk aged several million years is enough for the excitation process. 
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We also investigate the affection of disk mass on the exciting process. Figure [7] 
gives the scanning results like figure IHt for M^sk = 20Mj and M^k = lOOMj. Other 
parameters keep the same as in table 1 for simplicity. For the smaller disk mass (FigJTki), 
the regions of i2,o,min move toward upper and right relative to the same ones in figure [6^, 
which causes the zone of lower /2,o,min smaller. And for the bigger disk mass (FigJTb), the 
-^2,o,min = 5° ~ 7° range extends to the less massive M 2 -region as compared to figure and 
the /2,o,min = 0° ~ 5° range appears in the upper. So a more massive disk is in favor of the 
pumping to some extent. 

All the above discussions have set e2,o = 0.001 for concentrating on VSR more 
conveniently. However, m2 is more likely on an eccentric orbit since the planetary scattering 
have prompted a non-zero inclination. Figure [8] shows the same N-body simulations 
scanning as those in figure |3b,c with a higher e2,o- The remarkable difference in the higher 
e2,o situation is that the critical value of ^2,0 f° r pumping gets smaller. The cases with 
-fi.max > 90° and ei jmax > 0.99 appear even when 7 2 ,o = 0°, which might be due to the strong 
couplings between ESR and VSR when €2 is large. However, the effect is not obvious when 
e 2 , < 0.2. 

In Fig. 8, we also notice that at small e2,o (Fig.8a), planetary mean motion resonances 
(4:1, 5:1,6:1) cause an increase of -Zi jmax and ei jmax . When e2,o becomes larger, the regions 
outside 12AU in Fig. 8b and 8AU in Fig. 8c are full of the cases with Ji imax > 130° and 
ei,max > 0.99, this is due to that, as the apohelion of the inner planet is comparable to the 
perihelion of the outer planet, then planetary scattering dominates. We stop the simulation 
as long as any planet crosses the inner edge of the disk. 
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6. Systems with more than two planets 

With the help of the VSR, a mutual inclination between planetary orbits much smaller 
than the Kozai critical value can induce the pumping of the inner planet's eccentricity 
eventually. However, the occurrence of VSR constrains the inner orbit to a narrow trigger 
range, and the opportunity is small that two adjacent planets happen to be in the VSR 
configuration. Actually, the above pumping mechanism can be extended to multiple 
planetary systems so that a wider trigger range can be achieved. Here we give two examples 
to show that the mechanism can also occur between two nonadjacent planets, as well as 
inspire a chain reaction among more than two planets. In the left case of Figure the 
innermost and outermost planets were right in a configuration to be excited in a two-planet 
situation, and after another planet is added between them, the excitation still turns up. 
The right case in Figure |9] exhibits a chain reaction. The middle planet is located right 
in the VSR scope of the outermost planet, so its inclination is pumped at first, which 
directly leads to the increase of the mutual inclination of the inner two planets. At last, 
the innermost planet is excited by the VSR with the middle planet. So the influence of 
excitation of the outer planet can be spread to a more inward scope by a chain reaction. 
We do not explore the specific conditions or detailed influence of these more complicated 
operations of the mechanism, and leave them to future works. 

7. Conclusions and discussions 

In this paper, we proposed a mechanism to excite the eccentricities and inclinations 
of planets with a residual gas disk outside the planets. The excitation was the results 
of a coupling of secular resonance and Kozai effect. After several giant planets formed, 
the inner disk was assumed to have been swept out by gas giants during their accretion, 
and the outer part of gas disk would coexist with planets as long as million years. If the 
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outermost planet has a moderately inclined orbit relative to the disk midplane, vertical 
secular resonance would happen between the planets. Then the mutual inclination between 
two planets increases. Once it reaches up to ~ 40°, the Kozai effect between the planets 
would be induced, which can further pump the inner planet's eccentricity and inclination to 
high values (Fig. [1]). So this kind of mechanism is probably one of origins of hot-Jupiters 
on misaligned even retrograde orbits. 

To describe the evolution of inclinations and longitude of ascending nodes, we derived 
the evolution Equations, which are closed with the assumption of = 0, and are suitable 
for arbitrary inclinations. They are used to find out the locations and minimum initial 
inclinations for the occurrence of vertical secular resonance and Kozai resonance (Fig. [3j 
[6]). The elements here are relative to the disk midplane, and the formulas are different with 
those of the elements relative to the invariable plane of the two orbits. So they can be 
utilized to situations with the elements relative to any invariable plane. 

From the evolution equations, we showed that, with a mass of O.O5M residual gas disk 
located from 50au to lOOOau, a Jupiter-mass planet will be pumped by an outer gas giant 
with 5mj mass and 10° relative to the disk midplane at least, located out of 25au. And it 
could be flipped into a retrograde orbit by an outer gas giant with 10m j mass and an initial 
inclination of 20°. Such a mechanism can be also effective for a system with more than two 
planets, and the critical angles required might be more flexible with the presence of more 
planets. 

We used a simple disk model in order to compare with the results of the evolution 
equations and fully discuss the effect of planetary parameters. Also the limitation is that 
the disk mass have to be much bigger than the total mass of planets for satisfying the 
angular-momentum-advantage assumption(section 2), which restricts the full discussion to 
the disc parameters. To verify that the pumping process is irrelative to disk model, we 
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made the same simulations using a different disk model (INagasawa et al.ll2000t IZhao et al 



201ll ). Then we found the pumping still exists with similar structures and even locations of 



contour lines in Figure [3J 



The mechanism revealed above has some resemblanc e with that in binary systems. In a 



Takeda et al. 



(J2008|) divided three 



binary system with two planets orbiting one of the stars, 
distinct dynamical classes according to differential nodal precessions of the two planets. The 
mechanism illustrated in our paper is similar to the so-called "weakly coupled systems", 
with the same peculiarity that the planetary mutual inclination is excited by the secular 
resonance between the planets. That is actually a transitional case between "decoupled 
systems" and "dynamically rigid systems" . We illustrate this from the three cases in Figure 
[TUl where the semi-axis of the inner planet is the only varying parameter. In the left case, 
the planetary secular interaction is very weak and suppressed by the perturbation from 
the disk, so the secular nodal precession of the inner planet is much slower than that of 
the outer planet. In the right case, the mutual effects between the planets become so 
strong that their nodal precesses coupled, and the maximum of their mutual inclination is 
roughly the sum of ii,o and 1%$. The middle case is an exciting one, which occurs when 
planetary interaction is big enough that the secular nodal precessions of the two planets are 
approaching but not too big that the planets are coupled. 



In this respect, a protoplanetary disk has a comp arable effect wit h a ste llar companion. 

Wu fc Lithwickl ( 2011 ). They pointed 



Actually, this kind of analogy has been mentioned in 
out that the place of a planet in secular interactions could be replaced by a mass wire 
made by spreading the planet along its orbit. Since protoplanetary disks are universal in 
single-star systems, our exciting mechanism induced by secular resonance would not be 
limited in binary systems but can be extended to single-star systems. 



Though in our mechanism, eccentricity pumping can occur with an initial mutual 
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inclination much smaller than the Kozai critical angle, it is still within a rather narrow 
range of disk and planet configurations for a Jupiter-mass planet can be flipped. The 
efficiency for the occurrence of this mechanism in different systems will be investigated 
in future works. Comparing to observations, the narrow range also implies that, firstly, 
there should be most of systems owing planets with moderate or low eccentricities and 
inclinations than the systems owing retrograde hot-Jupiters. Secondly, according to our 
additional simulations, the more massive the inner planet is, the higher the initial outer 
inclination is demanded to be, meanwhile, the more massive the outer planet needs to be. 
So we speculate that the proportion of misalignment in Earth-like or Neptune-like planets 
is probably larger than that in Jupiter-like planet. All these need to be verified by further 
statistics of simulations as well as observations. 

The authors thank the referee for good suggestions which greatly improved this paper. 
The work is supported by National Basic Research Program of China (2013CB834900), 
Natural Science Foundations of China (10833001, 10925313), and Fundamental Research 
Funds for the Central Universities. 

APPENDIX 



A. Evolution of the orbital elements due to planetary perturbation 



We apply 



egendre polynomials expansion and Runge-Lenz vector introduced in 



Mardling fc Linl (120021 1 to deduce the elements' evolution due to planetary interactions. The 



quadrupole contribution of the acceleration of the inner orbit produced by the third body is 

Gm,2 



fi 



■v 



-(3xR-r) 



(Al) 



-18 - 



and that of the outer planet from the inner one is 



2,P 



i? 4 m 01 



-(5a; 2 — r 2 )R — 3xr 



(A2) 



where r and R are position vectors of the inner and outer planet in Jacobi coordinates, 
x — r ■ R, m i2 = m + m 1 + m 2 , m 01 = m + mi, fJ, i = m G m 1 /m m . 

The relations between the rates of change of the inner orbital elements and those of 
Runge-Lenz vectors are given by 

dei . „ 
e • e, 



dei 
~dt 

duji dVti e „ 

at at ei 

(i/i — (sin uii e + cos uiq) ■ h 
dt hi 

dQi (cosuie — smuiq) ■ h 
dt hi sin 

de = 2(/- r)r- (r- r)f- (/• r)f 

(it Gm i 

— = r x f 
dt J 

h = rxr is the orbital angular momentum vector of the inner orbit, e is the 
Runge-Lenz vector and q = hxe. r = ai(cosi£i — ei)e + ai a/1 — ef sin E\ q, r 



where 



(A3) 
(A4) 

(A5) 

(A6) 

(A7) 
(A8) 



— aiUi sin E\/ (1 — ti cos Ei)e + ainiyl — e 2 cos£?i/ (1 — ei cos £?i)qr. For the outer orbit, r 
would be replaced by i?, and the correspond unit vector is (E, Q, H). 

We separately substitute the expressions (I All) and (IA2I) for / in (IA3|) - (1A6|) . and average 
first over the inner orbit then the outer orbit, and simplify the results as follow 



dii 

[ dt )p 



3Gm 2 a\ , 



3 2n-3/2 



(cos LOi e — sin u)\ q) + e 2 (4 cos Wi e + sin q) 



[h-E)E+{h- Q)Q 



(A9) 
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,dVti. 

-IT- 



3Gm 2 a\ 
4hi sinii a 2 



3 (l-^)- 3/2 



(sin u>i e + cos U\ q) + e\ (4 sin u>i e — cos u\ q) 



(h-E)E+ (h Q)Q 



(A10) 



~dt 



)v = 



4:h 2 m 01 a s 2 
+ {l-e\)(H-q)q 



(1 - e\)- z/2 {cosuj 2 E - sinu 2 Q) ■ \{l + Ae\){H- e)e 



(All) 



,dVt 2 



3Gfi 01 m 012 a\ 
4h 2 m,Qi sinix a 2 



(1 - e 2 2 y 3/2 (sin tu 2 E + cos tu 2 Q) ■ \(l + Ae{){H- e) 



+ (l-el)(H-q)q 



(A12) 



,dei 
'~dt 



lhm 2 a\ 
Am m a 2 



n ien /l- e 2(l- e 2)-V 2 (e • E)(q- E) + (e - Q)(q- Q) , (A13) 



,duoi 



,dn 

~dl 



j 3m 2 a\ 
4m 01 a 2 



l-e\{l-el)-^{A[(e.Ef + (e.Qf 



\{q-Ef + {q.Qf}-2). 



(A14) 



The coordinates of the Runge-Lenz vectors relative to an arbitrary inertial plane are 

( cos^i cosc^! — sin^x sinwx cosii ^ 

sin cos uj\ + cos Vt x sin oj\ cos i\ 
y siniisinwi J 

( — cosfii sinwi — sinfii coswi cosii \ 
sin Qi sin u>i + cos Qi cos u>i cos i\ 

sinixcos^i J 

( siniisinfii \ 

— sinii cos Q\ 
\ cosii ) 



h 



(A15) 



and for E, Q, H, the formulas are similar except for switching the subscripts from 1 to 2. 



Then we derived the final simplified expressions for the rates of change of h,I 2 , ^1,^2,^1 
and uj\ 

?>m 2 a\rii 



, dh 



'dt )p 



4m ia 2 



1 - e 2 ) l,2 {l-el) 3/2 cos/i cos/ 2 + sin/isin J 2 cos(fii - fi 2 ) 



1 



x { sin I 2 sin(r2i — f2 2 ) + -e 2 (3 + 5 cos 2oo\) sin I 2 sin(f2i — f2 2 ) 



+5 sin 2oj\ (cos Ii sin I 2 cos^x — Q 2 ) — sin ^ cos I 2 ) 



(A16) 



3m mia 2 n 2 
4moi fl 2 



1-e 2 .) 



2\-2 



sin/i sin(r2! — Q 2 ) 
1 



COS /1 COS I 2 

3 sin /1 sin(fii — Q 2 ) ( cos ii cos 7 2 



+ sin ii sin 7 2 cos(f2i — Vt 2 ) 
+ sin/i sin J 2 cos^ — f2 2 )) + 5 cos2w! sin^x — (^(sinix cos fx cos/ 2 
— (1 + cos 2 Ii) sin I 2 cos(fii — Q 2 )) + 5 sin2w 1 ( sin I\ cos/ 2 cos(f2i — Q 2 ) 



— cos/i sin/ 2 cos 2(fii — f2 2 )) 



= 3m 2 a?m _ e 2)-i/2 (1 - e 2)-3/2J ± sin/ j 2cos2(^ 1 -^ 2 )sin 2 / 2 
4m ia^sin/i 14 L 



(A17) 



-3 cos 2/ 2 — 1 



I 1 

- cos 2/x sin 2/ 2 cos(^i — Q 2 ) + -e 2 



COS /1 COS I 2 



+ sin fx sin I 2 cos(f2i — Q 2 ) (—3 + 5 cos 2cui) ( sin Ii cos 7 2 
— cos Ii sin J 2 cos(f2! — f2 2 )) + 5 sin I 2 sin 2wi sin^x — f2 2 ) 



(A18) 



6^2 

1 dt )p 



3m mia 2 n 2 2 2 J 1 . 

1 — 2 2 • — r( 1_e 2) <7smi 2 cos/ 2 
4m 2 )1 a 2 ! sml 2 [4 

+ i sin 2/i cos 2/ 2 cos(fii — Q 2 ) + ^e 2 



2 cos 2(fti - ft 2 ) sin 2 h - 3 cos 2ii - 1 
3 sin 2J 2 ( — cos 2 ix + sin 2 1\ cos 2 (f2x — f2 2 )) 



+3 sin 2/i cos 2/ 2 cos(f2i — Q 2 ) — 5 cos 2wi sin 2ix cos 2J 2 cos(f2x — Q 2 ) 
+5cos2wi sin2/ 2 (cos 2 iicos 2 (f2i - Q 2 ) - sin 2 (f2x - Q 2 ) ~ sin 2 /i) 
+ 10 sin 2^ sin(f2x — Q 2 )( sin I\ cos2/ 2 — cos Ii sin2J 2 cos(f2x — VL 2 )) 



(A19) 



.dei 
■~dt )p 



,dui 



w 

of 



8moia 2 



3 2\-3/2 



sin2co' 1 



sin 2 J 2 sin 2 (Qi — Q 2 ) — 2 cos 2a>i sin J 2 sin^! — Q 2 
cos Ji sin J 2 cos(fii — Q 2 



sin cos I 2 — cos ii sin I 2 cos(£7i — £l 2 )) 2 
sin Ii cos J 2 

(A20) 



,dVLi ?>m 2 a\ni 



1 _ e 2(l _ e 2)-3/2 I (l_ 5sin 2 Wi; 



4moia 2 

+ cos Ii cos/ 2 cos(f2i — Q 2 )) 2 + sin 2 (f2i — f2 2 )( cos 2 1\ + sin 2 I 2 ) — 1 
—5 sin 2cui sin(f2i — Q 2 ) sin J 2 sin I\ cos J 2 — cos I\ sin J 2 cos(f2i — Q 2 

+3sin 2 / 2 sin 2 (fi 1 - Q 2 ) - 1 1. 



sin /i sin J 2 



(A21) 



len — = 71, t he latter two formula turn to the quadrupole parts of (C9) and (C5) 



Naoz et al. 



(2011b) 



B. Evolution of the orbital elements due to disk gravity 



As in observation, disk mass is a commonly estimated parameter rather than the radial 
distribution exponential or mass density, we set disk mass as an independent parament and 
deduce the mass density from 



R 



out 



2nrdr = M disk , 



then obtain 



with 7) = R in /R out . 



(-a + 2)M disk 

2(1 - T]~ a+2 )TTR 2 



(Bl) 



(B2) 



out 



We used Lagrange's equations in 



Murray fc Dermott 



(119991 ) to deduce the rates of 



change of elements due to disk gravity. First, we expanded the gravity potential in r p /r to 
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the quadrupole, like iTerquem et al.l (120 10f ) 



-a + 2 GM t 



disk 



1 _ v -a+2 R 



out 



1—1] 



1-a 



■1 + 7)- 



-1-a 



1 — a 



1 + a 2i? 2 



out 



-1+2 sin2 #p 



(B3) 



The first term in the square brackets has no contribution to derivation, so only the second 
one is retained. Defining 

-a + 2 -l + f'^GMrf, 



K 



'disk 



-1-a 2i? Q 3 ut ' 



then substituted the expresses with true anomaly / for r p and 9 P , we got 



(B4) 



2\2 



$ = AT- 



q 2 (l-e : , 
1 + e cos/) 2 



1 ^ • 2 / n\ . O , 

sin (a; + /) sin 1 

2 2 



(B5) 



We substituted the a bove one into Lagrange's Equations (6.148)-(6.150) in 



Murray fc Dermott 



finally 



(119991 ). then averaged over true anomaly /, and got the evolutions 

,da s 



dt 



J disk 



0. 



I disk 



,de s 
"dt 



"dt )disk 



sin 2oo sm 1 , 



4n 
15Ke 2 



sin2u; sin 21, 



jdisk 



3K cos / 

4n/3 



(2 + 3e 2 -5e 2 cos2cu), 



du K ( 9 2 15 2 o 

(_) dlsk =_j-2--e + T e cos2a; + sin / 

where (3 = a/1 — e 2 . 



4+16 6 16 (2 + e)GOS2a; . 



(B6) 
(B7) 
(B8) 
(B9) 
(BIO) 



When % ~ and e = 0, the expressions can be simplified into 



da. 

{ di )disk 



,de s 
dt 



,dL 



(disk 



disk 



dt 
du 

[ dt Jdisk 



"dt 
3K cos / 
2n 
2K 

n 



disk 



(Bll) 

(B12) 
(B13) 
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Table 1: Initial condition for Figured! 



Planet 


Mass 


Semimajor Axis 


Eccentricity 


Inclination 




(Mj) 


(au) 




(°) 


mi 


1 


3 


0.001 


1 


m 2 


10 


30 


0.001 


30 


disk 


Mass 


-Rin 


-Rout 


a 




(Mj) 


(au) 


(au) 






50 


50 


1000 


1 



Note. - 



Other arguments are taken arbitrarily except for longitude of nodes of two planets equal. 
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t (Myr) t (Myr) 

Fig. 1. — Evolutions of two planets with (left panels) / without (right) an outside disk's 
gravity. The initial conditions are listed in Table 1. Black lines are for the inner planet, and 
red for the outer one. Green lines in the top panels indicate the mutual inclination of two 
planets. Dash lines in the two left-upper panels and lighter dots in the two left-lower panels 
are the results of the evolution equations (jHJ). 
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t (Myr) t (Myr) 



Fig. 2. — The same as Figure [U except a^o = 5.5au,a2,o = 35.5cm, J 2 ,o = 10°. 
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Fig. 3. — Contours of maximum of the mutual inclination between two planets itot,max( a )> 
maximum of the inclination of the inner planet /i, m ax(b), maximum of the eccentricity of the 
inner planet ei iinax (c), maximum of the eccentricity of the outer planet e 2im a X (d) from full N- 
body simulations during the evolution of lMyr. Every point has different initial inclination 
of the outer planet ^2,0 {y axis) and different initial semi-major axis of the inner planet a^o 
(as axis). The black lines in panel a,b and c indicate the results of the evolution equations 
OH]), which are integrated 1 million years or truncated after e± > 0.99. The black stars in the 
two upper panels are used to label the positions whose coordinates are contoured in Figure 
(i 
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100 



10 r 



1 r 



0.1 r 




a^AU) 



Fig. 4. — Precession timescales of argument of pericenter u and longitude of ascending node 
Q of two planets. Initial parameter is listed in Table 1, except for i2,o = 0, and a\ altering 
from to 15au. Secular resonance for e — u> would take place around 3.3au, the place 
T un = T u) 2 i an d secular resonance for / — Q around 2.7au, the place tq 1 = tq 2 . 
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_4 I i I i I i Ll i I i 

0.0 0.2 0.4 0.6 0.8 1.0 

t (Myr) 



Fig. 5. — The evolution of \og[(dx/dt) p /(dx/dt)di S k\ (% represents I\, Qi, ei, u>i, I2 and Q 2 ) 
with time for the case in Figure [U Red line is the boundary where (dx/dt) p = (dx/dt)di S k- 



-32 - 




= 90° 



20 



15 



10 



(b) 


7 


"/ 

□ 












' o 
















I 










o 










<o 








o 

Q 


o' 








c 


I 




../ 




./ 








/ 





I 2 , , min (deg) 

■ 



I 
I 



5° 
7° 

0° 
20° 
30° 
40° 



30 40 
Oh ( au ) 



10 



20 



30 
o M (ou) 



40 



Fig. 6. — With different a 2j o and m 2 , the left panel displays the minimum of initial inclina- 
tion of the outer planet (filled color contour) for cases in which 7 tot could reach 40° during 
evolution (y coordinations of the star in Figure Eh)- The solid line contour is made up of 
locations of the inner planet when Jtot.max = 40° happens with the smallest 7/2,0 ( x coordina- 
tions of the star in Figured). The right panel has similar meanings except for I\ reaching 
90° during evolution (the coordination of the star in Figure (3b)- Every a^o — I2.0 scanning 
involved is from the same condition as the black line contours in Figure [3J 
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Fig. 8. — The same as the results of N-body simulation in Figure |3^ (the left ones) and|3b 
(the right ones), except for (a) e2,o = 0.2. (b) e2,o = 0.35. (c) e2,o = 0.5. The red numbers 
in the panel (a2) mean the period ratios of two planets. 
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Fig. 9. — Two cases of evolution of semi-major axis, inclinations and eccentricities of three 
planets, which orbit the center star with a disk outside. The left plot has three planets 
with mi = lmj,m2 = lnij, m3 = 5mj,ai = lOau, = 20au, 03 = 40au, I\ = l°,l2 = 
1°, I3 = 20°, and three planets in the right plot are mi = O.lmj, m2 = lmj, m% = 5mj, a\ = 
lau, a 2 = 10au, a 3 = 40au, Ii = 1°,I2 = 1°,-^ = 30°, The disk parameters are the same as 
Table 1. 
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t (Myr) t(Myr) t(Myr) 

Fig. 10. — These are three cases representing three different kinds of evolution of the inner 
planet. The only different initial condition is the semi-major axis of the inner planet a 10 , 
which is 0.56au, 4.33au and 9.4au from left to right. Other parameters are the same, mi = 
Imj, m 2 = 5mj, mdisk = 50mj, a 2 = 20au, R{ n = 30au, R out = lOOOau, I\ = 1°, I 2 = 31°, e\ = 
e 2 = 0.001,fii = Q 2 . All arguments are arbitrary. Black lines are for the elements of the 
inner planet and red lines for these of the outer planet. 



